/* benchmark guru interfaces */
#include <Utilities/RandomGenerator.h>
#include <Configuration.h>
#include <mpi/communicator.h>
#include <mpi/collectives.h>
#include <Utilities/Timer.h>
#include <OhmmsPETE/OhmmsMatrix.h>
#include <Message/OpenMP.h>
#include <benchmark/transpose.h>
using namespace APPNAMESPACE;

int main(int argc, char** argv)
{
  mpi::environment env(argc,argv);
  mpi::communicator mycomm;
  OhmmsInfo ohmms("transpose",mycomm.rank(),0,1);

  typedef  double real_type;
  int niters=10;
  int sizeN=4;
  int sizeM=4;
  bool debug=false;

  int ic=0;
  while(ic<argc)
  {
    string a(argv[ic]);
    if(a == "opt_s")
      sizeN=sizeM=atoi(argv[++ic]);
    else if(a == "opt_n")
      sizeN=atoi(argv[++ic]);
    else if(a == "opt_m")
      sizeM=atoi(argv[++ic]);
    else if(a == "opt_i")
      niters=atoi(argv[++ic]);
    else if(a == "opt_debug")
      debug=true;
    ++ic;
  }

  typedef double value_type;
  Matrix<value_type>  A(sizeN,sizeM), trA(sizeM,sizeN);
  for(int i=0; i<A.size();++i) A(i)=Random();

  //const int teng=DUMMY_TRANSPOSER;
  //const int teng=SWAP_TRANSPOSER;
  const int teng=MKL_TRANSPOSER;
  //const int teng=ESSL_TRANSPOSER;
  if(debug)
  {
    app_log() << "Before transposes " << endl;
    app_log() << A << endl;
    //Transpose2D<value_type,teng>::apply(A,trA);
    Transpose2D<value_type,teng>::apply(A.data(),A.rows(),A.cols(),trA.data());

    bool fine=true;
    const double eps=numeric_limits<value_type>::epsilon()*10.;
    for(int i=0; i<A.rows(); ++i)
      for(int j=0; j<A.cols(); ++j)
        if(abs(A(i,j)-trA(j,i))>eps) { cerr << "A("<<i << "," << j << ") != trA " << endl; fine=false;}
    if(fine)
      app_log() << "Transpose is correct " << endl;
  }

  Timer t1,t2;
  vector<double> dt(2,0.0);
  for(int i=0; i<niters; ++i)
  {
    //add random numbers
    for(int i=0; i<A.size();++i) A(i)+=Random();

    t1.restart();
    Transpose2D<value_type,teng>::apply(A,trA);
    dt[0]+=t1.elapsed();
  }

  for(int i=0; i<niters; ++i)
  {
    //add random numbers
    for(int i=0; i<trA.size();++i) trA(i)+=Random();

    t2.restart();
    Transpose2D<value_type,teng>::apply(trA,A);
    dt[1]+=t2.elapsed();
  }

  mpi::reduce(mycomm,dt);

  for(int i=0; i<dt.size(); ++i) dt[i]/=static_cast<double>(mycomm.size());
  double op_counts=sizeN*sizeM*niters*sizeof(value_type)/1.e9;

  //app_log() << "TRANS rows, cols, time(a-aT), gygabytes/sec, time(aT-a), gygabytes/sec " << endl;
  app_log() << sizeN << " " << sizeM << " " 
    << dt[0]/static_cast<double>(niters) << " " << op_counts/dt[0] << " " 
    << dt[1]/static_cast<double>(niters) << " " << op_counts/dt[1] 
    << endl;

  return 0;
}

